#!/usr/bin/python
# GM_SimpleIntegrate.py

#import math
import numpy
from time import asctime
from os import linesep

# Must have G, c in appropriate units so rE is in pc! Note M is in Msun
# See notes pg. 13 date 2/22/11
Goverc2 = 4.786e-14


def simple(simpars, phypars, appearance, masstuple, vtuple):
    '''Acquire lensing rates, qualified only by source magnitude.'''
    # Intake 'n' stuff
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    (darkstars, hz, dmlist, rho0mlist) = masstuple
    (vrelappl, sigmaappl2, sigmaappb2) = vtuple
    # Distance elements
    Distance = dlocs*PtScaled
    dDistance = delds*PtScaled
    DL = numpy.array(Distance)
    dDL = numpy.array(dDistance)
    DLdDL = DL*dDL
    DS = numpy.array(Distance)
    dDS = numpy.array(dDistance)
    DS2dDS = DS*DS*dDS
    del Distance, dDistance
    # Angular element, incl. symmetry restoration via quadrupling
    deltals = dells*PtScalela
    deltabs = delbs*PtScaleba
    bs = blocs*PtScaleba
    cosbs = numpy.cos(bs)
    cosbdeltab = cosbs*deltabs
    dOmega = numpy.multiply(deltals[:,numpy.newaxis], cosbdeltab[numpy.newaxis,:])*4 # Symmetry!
    del deltals, deltabs, bs, cosbs, cosbdeltab
    # These histograms compile things already "naturally" separated
    MvHist = numpy.zeros(len(limitset), dtype=float) # Divide by dMv here
    MassHist = numpy.zeros(len(MvMass), dtype=float) # Divide by dM here
    DSHist = numpy.zeros(len(dlocs), dtype=float) # Divide by dDS here
    DLHist = numpy.zeros(len(dlocs), dtype=float) # Divide by dDL here
    # These histograms need additional complex masks to delimit them
#    vtappHist =  numpy.ones(110, dtype=float) #### Not supported yet
    # This holds the overall result
    lensrate = 0.0
    lensrates = numpy.zeros(len(SourceMag))
    SourceMagIndex = 0
    CurrentLimitingMag = SourceMag[SourceMagIndex]
    for i in range(0, len(appearance)):
        nS = appearance[i]
        Mv = (limitset[i,0] + limitset[i,1])*0.5
        dMv = limitset[i,1] - limitset[i,0]
        # Cut for source unlensed magnitude
        if Mv > (CurrentLimitingMag + 0.001):
            lensrates[SourceMagIndex] = lensrate
            # If you've reached the end, stop; otherwise, move to next stage
            if SourceMagIndex == (len(SourceMag) - 1):
                break
            else:
                SourceMagIndex = SourceMagIndex + 1
                CurrentLimitingMag = SourceMag[SourceMagIndex]
        icontrib = 0.0
        for j in range(0, len(darkstars)):
            nL = darkstars[j]
            M = MvMass[j,1]
            ijcontrib = 0.0
            for k in range(0, len(dlocs)):
                Ki = k
                Kf = len(dlocs) - k
                # 1d
                two_rE = 2*(4*Goverc2*M*DL[:Kf]*(DS[Ki:] - DL[:Kf])/DS[Ki:])**0.5
                DLdDL_DS2dDS = DLdDL[:Kf]*DS2dDS[Ki:]
                AllDistances = two_rE*DLdDL_DS2dDS
                # 3d
                RMSvtapp = ((vrelappl[Ki:,:,:] - vrelappl[:Kf,:,:])**2 + sigmaappl2[Ki:,:,:] + sigmaappl2[:Kf,:,:] + sigmaappb2[Ki:,:,:] + sigmaappb2[:Kf,:,:])**0.5
                VolumeElements = numpy.multiply(AllDistances[:,numpy.newaxis,numpy.newaxis], dOmega[numpy.newaxis,:,:])
                GammaK = nS[Ki:,:,:]*nL[:Kf,:,:]*RMSvtapp*VolumeElements
                # Sum over b, then l, to have D-dependence
                GammaKblsummed = (GammaK.sum(axis=2)).sum(axis=1)
                DSHist[Ki:] = DSHist[Ki:] + GammaKblsummed
                DLHist[:Kf] = DLHist[:Kf] + GammaKblsummed
                ratecontrib = GammaKblsummed.sum()
#                print i, j, k, "ratecontrib", ratecontrib, asctime()
                lensrate = lensrate + ratecontrib
                ijcontrib = ijcontrib + ratecontrib
            MassHist[j] = MassHist[j] + ijcontrib
            icontrib = icontrib + ijcontrib
#            print i, j, "ijcontrib", ijcontrib, asctime()
#            print linesep
        MvHist[i] = icontrib/dMv
        print "For a width of", dMv, "around a source magnitude of", Mv, "I found a lensing rate contribution of", icontrib
    print "For a limiting unlensed source magnitude of", str(SourceMag)+ ", a lensing rate of", lensrate, "was found."
    # Divide each histogram by the width of the bin (i.e., normalize).  This has already been done for MvHist.
    MassHist = MassHist/dmlist
    DSHist = DSHist/dDS
    DLHist = DLHist/dDL
    return (lensrate, MvHist, MassHist, DSHist, DLHist)
                
